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We investigate how the time evolution of different kinetic Ising models depends on the initial 
conditions of the dynamics. To this end we consider the simultaneous evolution of two identical 
systems subjected to the same thermal noise. We derive a master equation for the time evolution of 
a joint probability distribution of the two systems. This equation is then solved within an effective- 
field approach. By analyzing the fixed points of the master equation and their stability we identify 
regular and chaotic phases. 
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The question to what extent the time evolution of a 
physical system depends on its initial conditions is one 
of the central questions in nonlinear dynamics that have 
lead to the discovery of chaotic behavior Q. In more 
recent years analogous concepts have been applied to 
the stochastic time evolution of interacting systems with 
a macroscopic number of degrees of freedom. Among 
the simplest of such many-body systems are kinetic Ising 
models where the above question has been investigated 
by means of so-called "damage spreading" simulations 
In these Monte-Carlo simulations two identical 
systems with different initial conditions are subjected to 
the same thermal noise, i.e. the same random numbers 
are used in the Monte-Carlo procedure. The differences 
in the microscopic configurations of the two systems are 
then used to characterize the dynamic stability. 

Later the name " damage spreading" has also been ap- 
plied to a different though related type of investigations 
in which the two systems are not identical but differ in 
that one or several spins in one of the copies are perma- 
nently fixed in one direction. Therefore the equilibrium 
properties of the two systems are different and the micro- 
scopic differences between the two copies can be related 
to certain thermodynamic quantities Note that in 

this type of simulations the use of identical noise (i.e. 
random numbers) for the two systems is not essential but 
only a convenient method to reduce the statistical error. 
Whereas this second type of damage spreading is well 
understood and established as a method to numerically 
calculate equilibrium correlation functions, much less is 
known about the original problem of dynamic stability. 
In particular, the transition between regular and chaotic 
behavior (called the " spreading transition" ) is not under- 
stood and it is unknown under which conditions it coin- 
cides with the equilibrium phase transition of the Ising 
model. (In numerical simulations it seems to depend on 
the dynamical algorithm whether the two transitions co- 
incide or not. Glauber dynamics usually gives a spread- 
ing temperature which is slightly lower than the equilib- 
rium critical temperature whereas the two are identical 
for heat-bath dynamics |q|.) 



In this Letter we therefore concentrate on the origi- 
nal question of the stability of the stochastic dynamics 
in kinetic Ising models. To this end we investigate the 
time evolution of two identical systems with different ini- 
tial conditions which are subjected to the same thermal 
noise. We derive a master equation for the joint proba- 
bility distribution of the two systems and solve it within 
an effective-field approach. By analyzing the fixed points 
of this equation we identify regular and chaotic phases. 
We find that the location of these phases in the phase 
diagram is sensitive to the choice of the dynamic algo- 
rithm. In particular, Glauber dynamics and heat-bath 
dynamics give very different dynamical phase diagrams. 

We consider two identical kinetic Ising models with N 
sites, described by the Hamiltonian 



- - - ^ JijSiSj - h^Si 



(1) 



where Si is an Ising variable with the values ±1. The 
dynamics is given by one of the two following stochastic 
maps which describe Glauber dynamics 



S',(i+ 1) = sign\v[h,{t)] 
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and heat-bath dynamics 



S.dt + 1) = sign {v[h,{t)] - Ut)} 



with 



vih) 



-h/T 



(2a) 



(2b) 



(3) 



Here hi{t) = J^j JijSj{t) + h is the local magnetic field 
at site i and (discretized) time t, ^i{t) G [0, 1) is a ran- 
dom number which is identical for both systems, and T 
denotes the temperature. Note that Glauber- and heat- 
bath algorithm differ only in the way how the random 
numbers arc used to update the configuration. The tran- 
sition probabilities v are identical for both algorithms. 

In order to describe the simultaneous time evolution 
of two systems H^^'> and ij'^^ with Ising spins 5*1^^ and 
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S'P'' we define a variable Vi{t) with the values ly = ++ 
for ^(1) = ^ 1^ fo^ ^(1) ^ _^(2) ^ 1^ _^ 

for -S-ti) = S'(2) = 1 and for S^^'> = S^^) = -1 

which describes the state of a spin pair (S'^^\ S'*-^^). Since 
we are interested in the time evolution not for a single 
sequence of ^^(i), but in ^-averaged quantities we consider 
a whole ensemble of system pairs {H'-^\H^^'> ) and define 
a probability distribution 



P{vi,...,VN,t) 



(4) 



where (•) denotes the ensemble average. The time evo- 
lution of P{vi, ...VN,t) for a single-spin dynamical algo- 
rithm as, e.g., Glauber or heat-bath dynamics is given by 
a master equation 

d , , 

— P(l/1, ...,l'N,t) = 

N 

^ P{vi,...,Vi,...,l'N,t)w{Vi—>^i) 
2—1 ^i^Vi 

N 

+ ^ ^ P(vi,...,^i^,...,VN,t)w{^i^->Ui). (5) 

Here w{iii ^ i^i) is the transition probability of the spin- 
pair (5'j-^\ S"-^^) from state ^ to v. It is a function of the 
local magnetic fields and h^^ and can be calculated 
from (2a) or (2b) for Glauber and heat-bath dynamics, 
respectively. 

A complete solution of the master equation (5) is, of 
course, out of question. Therefore one has to resort to ap- 
proximation methods, the most obvious being mean-field 
approximations. A natural way to construct a mean- 
field theory is usually to take the range of the interaction 
Jij to infinity. However, a mean-field theory constructed 
this way does not show the chaotic behavior found in 
the Glauber Ising model at high temperatures. A more 
detailed analysis shows that the absence of any fiuctu- 
ations in the infinite-range model is responsible for this 
discrepancy since the fiuctuations are essential for the 
chaotic behavior 

We therefore develop a slightly more sophisticated 
effective-field approximation that retains the fluctua- 
tions, though in a quite simplistic manner. The central 
idea is to treat the fiuctuations at different sites as sta- 
tistically independent. This amounts to approximating 
the probability distribution P(i/i, ...j/at, t) by a product 
of identical single-site distributions Py, 
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Piiyi,...,iyN,t) ^Y[P^^{t). 



(6) 



Using this, the master equation (5) reduces to an equa- 
tion of motion for the single-site distribution P^, 



j^P.^Y.^~P.W{,y^ fi)+P,,W{fi^iy)], (7) 



where 



(8) 



is the transition probability averaged over the states Vi 
of all sites according to the distribution P^. Note that 
the average magnetizations m^^^ , m^^^ of the two systems 
and the Hamming distance (also called the damage) 



i=l 



(9) 



which measures the distance between the two systems in 
phase space can be easily expressed in terms of P, 



m(i) = P++ + P+_ - P_+ - P__, 
m(2) = P_^+ - P+_ + P_+ - P__, 

i^ = P+-+P-+. 



(10a) 
(10b) 
(10c) 



So far the considerations have been rather general, to 
be specific we will now concentrate on a two-dimensional 
system on a hexagonal lattice with a nearest neighbor 
interaction of strength J. The external magnetic field h 
is set to zero. To solve the master equation (7) for the 
single-site distribution P we first calculate the transition 
probabilities v) between the states of a spin pair 

from one of the stochastic maps (2a) or (2b) and then 
average these probabilities over the states of the three 
neighboring sites of a certain reference site with respect 
to the yet unknown distribution P. This yields the trans- 
fer rates VF(/i v) which enter (7). The calculations 
involved are quite tedious but straight forward, they will 
be presented in some detail elsewhere 0. 

The resulting system of non-linear equations for the 

variables P++, Ph , P h> ^-nd P-- can first be used 

to calculate the thermodynamics. As expected, Glauber 
and heat-bath dynamics give the same results. In par- 
ticular, there is a ferromagnetic phase transition at a 
temperature determined by 



, 3J J 4 

tanh h tanh — = — , 

Tc Tc 3' 



(11) 



which gives Tc/J « 2.11. In the ferromagnetic phase the 
magnetization is given by 



I (tanh 3 J/T -I- tanh J/P) - 1 
I tanh J/P- i tanh 3 J/T 



(12) 



We now discuss the time evolution of the Hamming 
distance D between the the two systems which charac- 
terizes the stability of the dynamics. In contrast to the 
thermodynamics Glauber and heat-bath algorithms give 
very different results for the Hamming distance. We first 
consider the Glauber case. 
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The equation of motion of the Hamming distance can 
easily be derived from (7) and (10c). In the paramagnetic 
phase we obtain after some algebra 



dt 2^ 



2i?3)tanh— . 



(13) 



This equation has three stationary solutions, i.e. fixed 
points, D*, viz. £)* = which corresponds to the two 
systems being identical, = 1 where S^'^^ = —S^'^^ for 
all sites, and Dg — 1/2 which corresponds to completely 
uncorrelated configurations. To investigate the stability 
of these fixed points we linearize (13) h\ d = D — D* . 
The linearized equation has a solution d cx e~^* with 
Ai = A2 = itanh 3 J/T and A3 = -ifanhSJ/T. Conse- 
quently, the only stable fixed point is — 1/2. Thus, in 
the paramagnetic phase the Glauber dynamics is chaotic, 
since two systems, starting close together in phase space 
{D small initially) will become separated exponentially 
fast with a Lyapunov exponent Ai, eventually reaching a 
stationary state with an asymptotic Hamming distance 
D — 1/2. Note, that the Lyapunov exponent Ai goes to 
zero for T — > CX). Therefore the time it takes the systems 
to reach the stationary state diverges with T — > cxd, as 
has also been found in simulations 

We now turn to the ferromagnetic phase. In order to 
find the fixed points of the master equation (7) we can 
set the magnetizations of both systems to their equilib- 
rium values (12) from the outset. In doing so we exclude, 
however, all phenomena connected with the behavior af- 
ter a quench from high temperatures to temperatures 
below Tc. These phenomena require an investigation of 
the early time behavior and will be analyzed elsewhere 
0. For m^^^ = m(^) = m the equation of motion for the 
Hamming distance reads 



-D = -{D-2,D-' + 2D^) tanh — 
dt 2^ ^ T 



J 

T 



n? 2D tanh — - tanh 



tanh 



3J 



(14) 



This equation has two fixed points D* in the interval 
[0,1]. The first, r)\ — exists for all temperatures. The 
second fixed point D\ with < D3 < i exists only for 
T > Ts where the spreading temperature is deter- 
mined by 



3m tanh — = tanh — — . 



(15) 



This gives « 1.74 « 0.82Tc. The stability analy- 
sis shows that D\ = is stable for T < Ts and un- 
stable for T > Ts with a Lyapunov exponent Ai = 
i tanh 3 J/T - fm^ tanh J/T. The fixed point D* which 
exists only for T > is always stable. Consequently, 
we find that the Glauber dynamics is regular with the 
asymptotic Hamming distance being zero for temper- 
atures smaller than the spreading temperature Tg but 
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FIG. 1. Magnetization m, asymptotic Hamming distance 
D and Lyapunov exponent Ai as functions of temperature for 
the Glauber Ising model. Below Tc the curve for D has two 
branches corresponding to the two systems being in the same 
or in different free energy valleys. 



chaotic for T > Tg. Close to the spreading tempera- 
ture the asymptotic Hamming distance increases linearly 
with T — Ts which corresponds to the spreading transi- 
tion being of 2nd order. In contrast to the paramagnetic 
phase, where the two systems become eventually com- 
pletely uncorrelated, for Ts < T < Tc the asymptotic 
Hamming distance D is always smaller than 1/2 so that 
the two systems remain partially correlated. Directly at 
the spreading point the term linear in D in (14) vanishes. 
For small Hamming distances the equation of motion now 
reads dD/dt oc which gives a power-law behavior 

D{t) oc t^i. 



Analogously, for m 



(1) 



j(2) — ffjQ find two fixed 



points, D2 — 1 which exists for all temperatures and 



with ^ < Dl 



< 1 which exists for T > Tg only. is 



stable for temperatures T <Ts and unstable for T > Ts 
whereas Dl is always stable if it exists. The results for 
damage spreading in the Glauber Ising model within our 
effective-field approximation are summarized in Fig. 1. 

We now investigate the time evolution of the damage 
for the Ising model with heat-bath dynamics (2b). After 
calculating the averaged transition rates W{fi ^ ly) 0| 
and inserting them into (7), we obtain the equation of 
motion for the Hamming distance D. In the paramag- 
netic phase it reads 



dt 



-D 



3D 



tanh 



tanh 



3D^ r , 3J 

tanh — 

T 



3J 



, J 

tanli — 
T 



J_ 

T 



4/3 



(16) 



D^ \ 3J 

tanh — 

T 



3 tanh — 
T 



This equation has only a single fixed point in the physical 
interval [0,1], viz. D\ = [0. It is stable everywhere in 
the paramagnetic phase. Consequently, the asymptotic 
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FIG. 2. Magnetization m and Lyapunov exponent Ai as 
functions of temperature for the heat-bath Ising model. 



Hamming distance is zero for all initial conditions, and 
the heat-bath Ising model does not show chaotic behavior 
for T > Tc- The Lyapunov exponent is given by Ai = 
I tanh ^ + J tanh ^ — 1 < 0. The Lyapunov exponent 
goes to zero for T — + Tg, thus at the critical temperature 
we again find for small Hamming distances dD/dt oc —D^ 
which gives D{t) (xt~^. 

In the ferromagnetic phase for m^^^ = m^^-* = m the 
equation of motion is given by 



dt 



(17) 



3D 

T 

3D2 



(1 -|-m^)tanh 



3J 



(1 -3m2) tanh - -4/3 



, 3J 
tanh — 
T 



tanh — 
T 
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, 3./ 
tanh — 
T 



o 1 J 

6 tanh — 
T 



Here we also obtain only one fixed point Dl = which 
is stable for all temperatures. The Lyapunov exponent is 
given by Ai = Kl+rn^) tanh ^-h|(l-3rn2) tanh ^-1 < 
0. Thus, the behavior is not chaotic and the asymptotic 
Hamming distance is D = 0. Analogously, in the fer- 
romagnetic phase for to^^^ — — m^^^ = m we obtain a 
single stable fixed point I?2 — The results for dam- 
age spreading in the heat-bath Ising model within our 
effective-field approximation are summarized in Fig. 2. 

In conclusion, we studied the simultaneous time evo- 
lution of two kinetic Ising models subjected to the same 
thermal noise by means of an effective field theory. For 
the heath-bath dynamics we found that two only slightly 
different equilibrium configurations stay close together 
in phase space for all times in both the paramagnetic 
and the ferromagnetic phase, i.e. an equilibrated heat- 
bath Ising model does not show chaotic behavior. For 
the Glauber dynamics we found a richer behavior. For 
all temperatures smaller than a spreading temperature 
Ts the two equilibrium configurations stay together for 



all times. For T > T^, however, their distance increases 
exponentially which corresponds to chaotic behavior. In 
agreement with numerical simulations for Glauber dy- 
namics the spreading temperature T, is not identical to 
the equilibrium critical temperature but slightly smaller. 

As with any mean-field theory we have, of course, to 
discuss in which parameter region it correctly describes 
the physics of our system. Since we treated the fluc- 
tuations in a very simplistic way, viz. treating fluctua- 
tions at different sites as independent, our effective- field 
theory will be reliable if the fluctuations are small, i.e. 
away from the critical point. Therefore our theory cor- 
rectly describes the high- and low-temperature behavior 
whereas it might misrepresent the details close to the 
critical point. In particular, the questions under which 
conditions the spreading temperature coincides with the 
equilibrium critical temperature and whether the spread- 
ing transition is of 1st or 2nd order cannot be considered 
solved. Further open questions are connected with the 
influence of external magnetic fields, long-range interac- 
tions and disorder. Some investigations along these lines 
are in progress. 
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